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1. Introduction 



Models of interacting individuals can be understood as many-body systems of statistical 
mechanics, and tools developed originally in the context of physics may be employed to 
address their dynamics and stationary states. This approach has fruitfully been apphed 
to a variety of agent-based models inspired by economics and game theory, see e.g. the 
recent textbooks [H [21 [3]. Attention here focuses on the interplay of co-operation and 
competition between interacting agents, and on the efficiency of their use of external 
information and resources. Statistical mechanics here offers a variety of valuable tools to 
study the global co-operative behaviour of such systems, and to understand their phase 
structure. In particular disordered systems theory [1] allows one to address interacting 
agent-models in which interaction matrices are drawn from random ensembles, and to 
compute typical average quantities for such models. Real-world systems are of course 
not random, but highly correlated. The aim of statistical mechanics approaches is hence 
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often not to study specific instances, but rather the general properties of classes of models 
as a function of the parameters characterising the distribution from which couplings are 
drawn. One may ask for example whether quantities such as connectivity, homogeneity 
or the strength of interaction affect the stability of a given model system. Taking an 
ensemble average here in a sense corresponds to studying all possible realisations of a 
given model at the same time, and hence to making statements about effects of model 
parameters in general, as opposed to analyses of specific real-world instances. 

This approach has been used to study e.g. the effects of self-interaction and 
memory in models of financial trading [H [2] or to examine how co-operation pressure, 
order of interactions impacts on the stability and trajectories of replicator systems of 
evolutionary game theory O [6|, [7] . In the context of population dynamics models with 
random interactions were first addressed by May in [8]. 

In this paper we study a model of a simple food-web composed of species and 
resources, originally proposed in a more basic form in [9]. The level of the resource 
consumption by species and its relationship with the stability of the ecosystem and 
the species richness is one of the main issues in ecology [10]. In [9], interaction between 
species is not through direct interaction (e.g. via prey-predator relations) but exclusively 
through the use and dependence on resources. If for example species A consumes a 
resource which B feeds upon as well, then this introduces a negative and symmetric 
interaction between A and B. The strength of negative interaction between A and B is 
hence regulated by the overlap in their dependence of resources. 

Due to the symmetry of interactions the discussion of p] focuses on a static analysis 
of this model eco-system. Here we choose a complementary dynamical approach, which 
allows us to address a broader class of interaction modes. Static studies necessarily 
rely on the existence of a Lyapunov function, extremised by the trajectories of the 
ecosystem, and are hence limited to systems with symmetric interactions. In the case 
of an ecosystem this is an obvious drawback, as competitive interaction of e.g. prey- 
predator pairs can not appropriately be addressed. A direct study of the dynamical 
equations allows us to extend the analysis to cases of asymmetric interaction matrices, 
and in particular to discuss the effects of anti-correlation on the behaviour of the system 
[TT| IT2] . Asymmetric interaction come in two ways in the present eco-system. Firstly, 
we introduce direct interaction between species, in addition to the indirect interaction 
through the use and dependence on resources. Secondly, we study the effects of possibly 
asymmetric dilution of the network of interacting species. 

The aim of our work is here twofold. Firstly, the study of the present model extends 
the statistical mechanics analysis of existing replicator models [5l El [3,1111 [12], relates 
to studies of Minority Games [Il[2]. Complex phase behaviour and different patterns of 
ergodicity breaking and instabilities have been identified in such models, with similarities 
as well as differences between rephcator-type models and other systems. One purpose of 
the present work is thus to contribute to the classification of such models according to 
the different types of phase transitions they exhibit, and to identify possible universal 
features. 
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Secondly, the model system here has a clear-cut ecological interpretation, even 
though the model may be criticised for not accurately capturing many features of real- 
world eco-networks. While our approach is a dynamical one, and ultimately results in a 
stochastic process for a single 'effective' species, all disorder in our model is quenched, 

1. e. the interaction web and coupling strengths are fixed at the beginning, and replicator 
dynamics are then considered on this fixed network. This approach on the one hand 
makes the model analytically tractable and allows one to reduce the description to a set 
of a few non-linear equations describing the relevant order parameters. One the other 
hand it constitutes a considerable restriction with respect to real-world eco-system, in 
which the web of interaction is of course not fixed, but subject to dynamical evolution 
itself, requiring the study of the dynamics of the network itself in combination with 
population dynamics on the network. Such evolving food-web models have for example 
been presented in [131 [H [IS [Ml [HI [IE] . Related work is also found in [IHl EQl El]. 
Results here rely mostly on numerical simulations (see however [T3] for descriptions in 
a Master equation formalism) and the food-webs resulting from these models have been 
compared to real-world data with respect to quantities such as the number of trophic 
levels, their relative populations and the typical connectivity of species. These models, 
some of which combine initial Gaussian random score matrices with evolving species 
networks, clarified the necessary conditions of types of functional responses and dietary 
choices (specialist/generalist) for producing realistic webs, whose structure agreed with 
empirical data. 

From the technical point of view it is interesting to note that recent stochastic 
models of complex food-webs [TU [15] and the 'neutral' model [22] effectively reduce 
multispecies stochastic process to a 'one species' process of a representative species which 
is subject to a 'mean-field' interaction with the remaining system, and that these models 
derive reasonable species abundance distributions in good agreement with real data. In 
a similar fashion our approach reduces the dynamics of species randomly coupled via 
quenched interaction to a 'one species' effective process as well. This mapping is fully 
exact in the thermodynamic limit in the statistical sense. Apart from providing a 
starting point for more realistic modifications of the present model, our analysis can 
hence, to a certain degree, be seen as complementary to the approach of [Ml [15]. 

The paper is organised as follows: we will first define the model, and then briefly 
discuss the statistical mechanics analysis based on a path-integral approach. We 
then turn to a stability analysis, and then discuss the effects of resource variation, 
direct interaction between species, co-operation pressure and dilution in the subsequent 
sections. We summarise our results in the conclusions section and point out lines for 
potential future research. 

2. Model Definitions 

The model describes an eco-system consisting of species, labelled by z = 1, . . . , 
and P = aN resources /x = 1, . . . , c^N. a is here a model parameter and is taken not 
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Figure 1. Illustration of the model: species compete for scarce resources while at the 
same time being subject to direct interaction e.g. through prey-predator relations. 



to scale with N, i.e. we assume a = 0{N^). The composition of the eco-system at 
time t is described by concentrations Xi{t) of species i = 1, . . . ,N, which evolve in time 
according to the following replicator equations [23] 



fi here denotes the fitness of species i at time t, and is frequency dependent. To 
be more precise fi is taken to depend on the composition of the ecosystem x(t) = 
{xi{t), . . . ,XN(t)) as well as on the abundance of resources Q(t) = {Q^(t), . . . ,Q^{t)). 

is a global 'field' variable, which is (up to a sign) typically chosen as the mean 
fitness in order to maintain the overall concentration of species. 

We will in the following assume that the fitness of species is composed of three 
contributions fi(t) = /j,s[x(t)] + /i,r[Q(t)] + fi^c{xi(t)). fi^s denotes a term describing 
direct species interaction, fi^r refers to interaction due to competition for resources. 
These two components of the model are illustrated in Fig. [H and can be understood 
similar to what is referred to as basal and intermediate species for example in [16]. a 
thus controls the relative number of basal species (resources) over intermediate species 
in our model. Finally /j c is an additional contribution describing an external so-called 
cooperation pressure, driving the eco-system to a state in which all species are present 
at equal concentration. We will in the following detail these three contributions to the 
fitness further. 

Following [6] we will choose the direct interaction between species to be 
characterised by a random couplings, i.e. 



Xiit) 



/,[x(t),Q(t)] + i.(t). 



(1) 



Xiit) 



N 




(2) 
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where the matrix elements {wij} are chosen from Gaussian ensembles according to the 
following distribution 

„, > iV / iV(^.g,-2^^t■,j^.,i + tt■|.) ^^ 



for any pair i < j. The diagonal elements are taken to vanish, wu = 0. Denoting the 
average over the random couplings by an overbar one thus has 

wij=0, wl = —, w-w-=T—. (4) 

r is a symmetry parameter and takes values F G [—1, 1]. For F = 1 the interaction 
between any pair of species i < j is fully symmetric, Wij = Wji. For F = Wij and 
Wji are uncorrected, and F = — 1 corresponds to a prey-predator relation, Wij = —Wji. 
Choosing intermediate values of F allows one to interpolate smoothly between these 
regimes. The ecologically most relevant setup presumably corresponds to negative values 
of F, describing competitive direct interaction between species, rather than co-operation. 

The second contribution fi^r to the fitness of species i describes its propensity 
to reproduce due to the presence or otherwise of resources. We here follow the lines 
of [9]. Let us assume that the amount by which species i G {1,...,A^} relies on 
resource /i G {1,...,P} is described by a coefficient ^f, with large signalling a 
strong dependence of i on /i. Then we will take /i,r[Q] to be of the form 

In turn a large abundance of i will then deplete the abundance of fi so that we write 

g'^(t) = g^[x(t)] = g{;-5^e^.(t). (6) 

j 

Qq here denotes the abundance of resource fj. in absence of species and the second term 
on the right-hand side corresponds to the consumption of resource n by the different 
species j = 1, . . . ,N. Recall that large indicates that species j consumes resource /x at 
a high rate, thus a large concentration Xj{t) (equivalently, a large number of individuals 
of species j) adds to the depletion of resource /i. The availability Q^(t) of resource /i 
thus becomes time-dependent, as the concentrations of species {xj{t)} evolve in time. 
In particular it appears interesting to ask the question whether or not the system is able 
to organise in a state which avoids over- and under-exploitation of resources, i.e. a state 
in which all Q^(t) remain close to zero asymptotically. We will address this question 
below. Following our earlier approach we take the coefficients {^f } to be drawn from a 
random distribution, specifically we choose them to be independent Gaussian variables, 
with mean q and unit variance, i.e. 



er = q, (er)^ - (O' = i- (7) 

According to the above remarks they describe the interaction between the species layer 
of the eco-system and the resource layer. While the following analysis focuses mostly 
on the case of Gaussian {^f} the generating functional theory below and computer 
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simulations show that only the first two moments of the {$,^} are relevant, so that more 
general distributions can be addressed as well with the methods used here. To complete 
the definition of fi^r, it remains to specify the {Qq}- Following [9] we write 



with {C^}^t=i,...,p independent standard Gaussian variables. The model parameter a thus 
controls the variability of resources. The scaling with P = aN of the {Qq} is chosen to 
guarantee a well defined thermodynamic limit, with which the theoretical analysis will 
eventually be concerned. 

Finally, we will study the effects of co-operation pressure on the eco-system. This 
variables acts to suppress the growth of individual species and is incorporated by a 
contribution 



to the fitness of species i [HI El] . In an ecological setting u takes mostly positive values 
denoting intra-species competition (but see also a comment on potential settings with 
negative u below). For u ^ oo the ecosystem is found in a state of perfect co-operation 
and maximal diversity (with all species surviving and having equal concentrations). As 
we will confirm later, a reduction of u leads to a smaller number of surviving species, 
and hence a reduced diversity. In order to obtain a complete overview of the phase 
behaviour of the model, we extend the analysis to negative values of u as well. 

The definition of the dynamics ([T]) is completed by stating the choice of z/(t) we 
will make in the following. In the analysis of the statics of the model it was found that 
only states with the normalisation A^~^ "^i^i = <^/(l contribute to the thermodynamics 
of the system |9J. Accordingly, we also restrict the dynamics to such configurations, and 
choose initial conditions and the subsequent Lagrange parameters {i^(t)}t>o such that 
the constraint 



is fulfilled at all times. This amounts to the choice uit) = —^j^ J^i^ii't) fil^it)^ Q(^)] • 
To conclude the presentation of the model let us briefly point out some of its obvious 
limitations. Firstly, due to the Gaussian choices of the {^f } and of the {Qq} negative 
values of these quantities might statistically occur (in the cases of the abundances {Qq} 
this is however suppressed in the thermodynamic limit due to the scaling with in ([8])). 
Secondly, the replicator dynamics ([I]) do not guarantee that all Q^{t) remain positive at 
all times. These drawbacks are consequences of the solubility of the model, as models 
with non-Gaussian disorder at the same scaling with A^ or additional constraints on the 
resource abundances are difficult to treat analytically. We would however like to note 
that with our choice of parameters (e.g. g = 1) most of the {^f } are indeed positive. The 
model is furthermore invariant under simultaneous shifts of the means of all {Qq-, } so 
that their averages can be chosen sufficiently high as to minimise the amount of negative 
couplings. 




(8) 




(9) 




(10) 
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3. Generating functional and effective species process 
3.1. Effective macroscopic theory and fixed point ansatz 

The model lends itself nicely to the study by the tools of disordered systems theory. For 
fully symmetric couplings F = 1 one identifies 

as a Lyapunov function, minimised by the replicator dynamics ([I]). Thus the stationary 
state of the model can in this case be obtained by purely static considerations based on 
replica theory. For general symmetry F no such Lyapunov function can be found, and 
the analysis needs to deal directly with the microscopic dynamics. The method of choice 
is here based on generating functionals, originally proposed in the context of random 
replicators in [6] , and recently used in [HI [T2] . The analysis focuses on the dynamic 
partition function 



Z[^] = ( (exp z / dtJ2^^m^{t) ) ) (12) 




where the average ((•)) extends over all trajectories of the system permitted by the 
equations of motion. 'J' is a source field introduced to generate dynamical correlation 
functions, and Z[^] is hence the Fourier transform of the probability measure on the 
space of paths generated by the replicator equations. Z[^] can then efficiently be 
averaged over the disorder, and evaluated by the method of steepest descents in the 
thermodynamic limit oo. We will not enter the detailed mathematics here, but will 
only report the final outcom^. One finds a description in terms of effective single-species 
trajectories, described by the following multiplicative Gaussian stochastic process 

x{t) = x{t)( [ dt'R{t,t')x{t') -ri{t) + u{t)). (13) 



^ Jto 

(to denotes the time at which the dynamics is started). The key components are the 
retarded interaction kernel 

R{t, t') = -2uS{t - t') - Tw'^Git, t') - a{l - G)-\t, t') (14) 

and the coloured Gaussian noise {ri{t)} which exhibits temporal correlations of the form 

{V{t)r]{t'))^ = w^C{t, t') + a [(I - G)-\aa^E + C)(I - G^Y^] (t, t'){lh) 

The matrices C and G in ((Hj) and (fT5l) are the correlation and response functions 
of the system, respectively, and are to be determined self-consistently as 

C(t,t') = (x(t)x(t'))., ^(^'0 = -(^) ' (16) 

I Imposing the above normalisation pop ensures that no super-extensive terms are found in the 
generating functional analysis and that the usual saddle-point integration can be carried out in the 
thermodynamic limit. 
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where (■)^ denotes an average over trajectories of the effective stochastic process ( IT^ . 
i.e. over reahsations of the noise {ri{t)}. E is the matrix with all entries equal to one, 
E{t,t') = 1 for all t,t'. 

The analysis proceeds by making a fixed point ansatz x(t) = x,ri{t) = r], u{t) = v 
in the effective process, leading to 

= (17) 

Furthermore we assume time-translation invariance of the response, i.e. G(t^ t') = 
G{t — t') and define the integrated response as 

X = j dtGit), (18) 

which we require to be finite for the further analysis. This restricts the ansatz to the 
ergodic regime of the system, i.e. to model parameters for which the assumed fixed- 
point is independent of initial conditions. The following self-consistent equations of the 
persistent order parameters {Q, x, J^} can then be derived along the lines of [HI [lH [T2] : 

(2u + w^Vx + -^^^ = [ Dz{A-z), (19) 



a 



— oo 



^2u + w^Tx + ^-) =[ Dz{A~zy, (20) 



A V 1 - X 



— oo 

-(2u + w''Tx + T^]x = I Dz. (21) 



oo 



Here Dz = ^=e ^ ^'^dz denotes the standard Gaussian measure, and one has A = w'^Q + 

V 2n 

a{aa^ + Q)/(l - x? and A = u/VX. We note that = f^^Dz = + erf (A/v^)) 
describes the fraction of surviving species. 

3.2. Key observables 

We will in the following study the behaviour of the system as a function of the different 
model parameters and in particular focus on the effects of the different components in 
the setup of the ecosystem. The above theory allows us to compute the behaviour of the 
model in the stable fixed-point regime exactly in the thermodynamic limit, and to carry 
out a linear stability analysis to identify the onset of instability as described below. 
Theoretical results will be compared to observations in computer experiments based 
on a numerical integration of the replicator equations ([T]). We here use the scheme 
of [7], effectively corresponding to a first order integrator with dynamical time step. 
In addition to the above mentioned fraction of surviving species 0, we will study the 

2 

diversity index D = j^q^, closely related to what is known as Simpson's diversity index 
in ecology [25]. Note that if the species concentrations were normalised to one the sum 
Ylii^f (i-6- the analogue oi Q = N^^Ylii^D would indicate the probability that two 
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randomly chosen individuals belong to the same species. We will also focus on the 
effectiveness of the use of resources. To this end one defines 

^=i^E((^"W)*)' (22) 

with (■)^ a time average in the stationary state. Note that {Q''{t)) = 0{N^/'^), so that 
H is of order one in the thermodynamic limit. H denotes the efficiency with which the 
species make use of the resources present in the system. If if = then {Q^{t))^ = 
for all /i, i.e. all resources are optimally exploited. If however H > 0, then the use of 
a fraction of resources (those with {Q^{t))^ ^ 0) is unbalanced. Our analytical theory 
allows us to compute H from the saddle-point equations, and one finds 

as in p]. 



3.3. Stability analysis and phase transitions 

The above ansatz of a stable ergodic regime breaks down, when either fixed points 
become numerous or suppressed in the thermodynamic limit. In the first case the system 
has a large number of (possibly marginally stable) attractors, and initial conditions 
determine which of these is realised. Hence ergodicity is broken. In the second case the 
system would not evolve into any fixed point at all at long times. 

The breakdown of the fixed point regime can be identified by means of linear 
stability analysis on the level of the effective process. Details of similar calculations 
can be found in [6l [11] . For the present model one finds that system runs into a unique 
stable fixed point if 

wY + " (TZ^ < (24) 

and that it becomes unstable when this condition is violated. 

Our above fixed-point ansatz also implies the assumption that the integrated 
response x be finite. A singularity in x would hence signal the breakdown of the 
ergodic theory and the onset of memory effects, in the sense that perturbations in the 
stationary state do not decay, but remain permanent p]. Simultaneously, a divergence 
of X necessarily implies H = (see eq. ( |23l) ). and hence a transition to a phase in which 
resources are optimally exploited. Since the right-hand-side (the fraction of surviving 
species) of Eq. ( |2T1) is bounded (0 G [0, 1]), we find that a divergence of x can occur 
only if M = and w = 0. Thus we expect no phase with optimal resource exploitation 
whenever co-operation pressure or direct species interaction are present. Finally we note 
that ( |2T1) implies = aif|x|-^c)oina model system with u = w = 0. Thus, the 
instability condition fl24|) is violated whenever x diverges. 
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4. Effects of resource variability 

We now first examine tlie effects of the variability of the resources. To this end, we 
set the strength of the direct species interaction w and the co-operation pressure u to 
zero in this section, and focus on the behaviour of the model as a function of a. This 
control parameter a measures the fluctuations of Qq (see Eq. ([H])), i.e. the degree to 
which the different resources fi = 1, . . . , C(N vary in their bare abundances Qq in the 
absence of species. For simplicity we keep g = 1 throughout this section. This system is 
the model studied in [9] by static methods. A phase transition was found, marked by a 
divergence of the static susceptibility in a replica symmetric ansatz. We here reproduce 
this transition from a dynamical calculation, and present the results of this section 
mainly for completeness and to set the scene for the subsequent parts of the paper. 

Fig. [2] shows the phase diagram obtained by solving the three equations fll9|20|2ip . 
The transition point ac = ac(o") is identified as the point where the integrated response 
X diverges. To obtain an interpretation of this transition in terms of the ergodicity 
properties of the system, we run two copies {x(t)} and {x'(t)} of the system with the 
same realisation of the disorder, but started from different random initial conditions 
and measure the distance (f = {N~'^ J^ii^ii'^) ~ ^'i{^))'^)t between two stationary states 
of the system. Thus if c?^ = initial conditions play no role, while for d? > {] the 
system is sensitive to the starting point. Although numerical measurements of d"^ can 
exhibit finite-size effects, simulations shown in Fig. [3] are consistent with an ergodic 
phase above Oc, and with a phase in which the system is sensitive to initial conditions 
below ac- In this second phase the system is still found to evolve into a fixed point, 
but stationary points of the dynamics become numerous, and which one of these is 
reached asymptotically is determined by initial conditions, similar to the behaviour of 
other replicator systems O [6l [7] . Fig. H] shows that this ergodic non-ergodic transition 
coincides with a transition between a resource-efficient phase at a < etc {H = 0) and an 
inefficient phase {H > 0) in the phase at a > ac- 

In Fig. [5] we report on the diversity of the eco-system as a function of the resource 
variability. One finds that the diversity of the ecosystem is large at a large number 
of resources per species, and that the ecosystem becomes less diverse as the number 
of resources is reduced. The figures also confirm that the behaviour of D = a'^/{q'^Q) 
is similar to the one of the fraction of surviving species (p, hence verifying the role of 
D as a measure of the diversity of the ecosystem. In the following sections we will 
hence focus on 0. As anticipated in the introduction, results do not depend on the 
specific shape of the distribution of the as only their first and second moments 

enter in the derivation of the effective dynamics. We have explicitly confirmed this in 
simulations, which show that measurements of H and D of systems in which the {^f } 
follow fiat, exponential, bimodal and power-law distributions with suitable first and 
second moments fall precisely on the lines obtained from the theory in Figs H] and [51 

The left panel of Fig. [5] furthermore confirms that = a at the transition with 
diverging susceptibility x- Similar transitions in static contexts can be identified by 
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Figure 2. Phase diagram for the model without direct species interaction and in 
the absence of co-operation pressm'e in the {a, a) plane {q — 1). The integrated 
response diverges at the phase transition line, and the system becomes sensitive to 
initial conditions in the imstablc phase. 




Figure 3. (Colour on-line) Reduced distance cP /Q versus a for the model without 
direct species interaction and co-operation pressure (w = O^u = 0, q = 1). Symbols 
show results from simulations for N > 200 species, run for > 10000 discretisation steps 
and averaged over at least 20 samples of random resource consumption {^,;}, vertical 
dashed lines mark the location of the phase transition as predicted by the theory. 



the divergence of the static susceptibihty in a rephca symmetric approach [9j. The 
occurrence of this transition has a geometrical interpretation similar to what is know in 
the context for example of Minority Games [23 EZl [H [2] ■ In the absence of co-operation 
pressure and direct species interaction, the fitness in expression ([1]) is of the form 
/i[Q] = J2fj.Q'^{'t)^i y a linear combination of the P A^-dimensional vectors 
= (^f, ■ ■ ■ ,'Civ)- "^^^ dynamics of the system thus can only wash out perturbations 
within the space spanned by the aN vectors = 1, . . . , c^N. Disregarding the 

(1 — (f))N extinct species, the underlying dynamical system has (pN effective degrees 
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Figure 4. (Colour on-line) H versus a for the model without direct species interaction 
{w = 0,q ^ 1). Curves are for a — 2,1,0 from top to bottom. The solid lines are from 
the theory, symbols from simulations {N = 300, 50 samples, run for 20000 steps). H 
vanishes below ac- 
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Figure 5. (Colour on-line) Diversity parameter D = a^/lQq^) (left) and fraction 
of surviving species (j) (right) versus a for the model without inter-species interaction 
{w — 0,q — 1). Curves are for cr = 2, 1, from bottom to top. The solid lines are from 
the theory, continued as dashed lines into the unstable phase in the left panel. Dashed 
line in right panel marks (j) = a. Symbols from simulations (parameters as in Fig. 3]). 



of freedom. Extinct species are typically stably extinct with respect to perturbations, 
see also [6]. The space of all potential external perturbations is hence 0A^-dimensional. 
Thus if a < some of those perturbations can not be removed by the dynamics, and 
ergodicity breaking occurs. 

The existence of a phase with H = at a < ac can be interpreted similarly. In the 
absence of direct interaction and co-operation pressure one has H and H coincide up to 
pre-factors, and the dynamics minimises this Lyapunov function. Attaining the absolute 
minimum H = implies {Q^{t))^ = for all fi via (!22l) . This constitutes a system of 
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aN constraints. With (j)N effective degrees of freedom available, these conditions can 
be met if a < but not above the transition point defined by (f){ac) = olc- 

5. Effects of direct species interaction 

Animals have not only the resource competition but also have the direct species 
interaction, like prey-predator, co-operation, competition and so on. In this section 
we study the effects of direct interaction between species, as controlled by the model 
parameter w. In order to focus on the impact of this model parameter we set m = 
throughout this section. We also limit the discussion to the case g = 1. 

In Fig. [6] we depict the phase behaviour of the system as direct species interaction 
is introduced. As predicted by the theory we find that the integrated response x is 
finite in case w > for all tested values of the model parameters, the transition lines in 
Fig. [6] hence mark an instability at which H remains positive. They are obtained from 
Eq. (!24l) . The figure demonstrates that the phase diagram is indeed affected by the 
direct interaction between species and by the symmetry of the couplings {wij\. The left 
panel shows that even a relatively moderate direct interaction of strength w = 0.1 can 
have a significant effect: symmetric (F = 1) and asymmetric (F = 0) interaction reduce 
the stable area while antisymmetric (F = —1) interaction expands the stable region, 
compared with the case without direct interaction [w = 0) shown in Fig. [2l This is 
confirmed in the right panel of Fig. [6j For symmetric and uncorrelated interaction 
ttc increases with increasing w so that direct interaction tends to make the system 
increasingly less stable. For negatively correlated interaction (F = —1) on the other 
hand, etc is a decreasing function of w, indicating that prey-predator-type interactions 
stabilise the ecosystem. One might speculate that for that reason, food-webs with this 
type of interaction may be more likely to be observed in nature than others. For F = — 1 
we also find ac approaches zero for large values of w indicating that there is no unstable 
region in the limit of w oo, which is consistent with marginally stable dynamics in 
the antisymmetric random replicator model without resource competition |28j . 

The left panel of Fig. [7| shows that the efficiency of resource exploitation is reduced 
as direct species-interaction is introduced, and is consistent with the predicted absence 
of a phase in which H = 0. The effect is stronger for correlated couplings than for 
negatively correlated ones. The effects of the direct couplings on the diversity of the 
eco-system is shown in the right panel of Fig. [71 One observes relatively little effect for 
the case of antisymmetric couplings, but a strong reduction of diversity as uncorrelated 
or positively correlated couplings are introduced. Crucially we here find that H and (j) 
are smooth functions of w as long as a > adw = 0). In particular no singularities are 
observed as w; — >■ 0. This is different in the case a = 0.2 < adw = 0) ~ 0.27, as shown 
in the insets of Fig. [71 Here i7 — > as w ^ and the integrated response diverges. 
Simulations at finite reveal non- monotonous behaviour of at w = 0"^. While we 
cannot fully rule out finite-size effects similar discontinuities of order parameters have 
been found in the context of so-called grand canonical Minority Games [H, [2]. The 
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Figure 6. (Colour on-line) Effect of direct species interactions. Left: Phase 
diagram in the (<t, a) plane {w — 0.1, q — 1, u = 0). The curves are obtained from 
Eq. (|24|) and are shown for F = 1,0, —1 from top to bottom. System is stable above 
the respective curves, and unstable below. Right: Phase diagram in the [w, a) plane, 
cr = 1, g = 1 and u = 0. Curves are for T = 1, 0, —1 from top to bottom. System is 
stable above the respective curves. 




Figure 7. Effect of direct species interactions: efficiency of resource exploitation 
H (left) and fraction of surviving species (j) (right) versus w {u — 0, a = q = a = 1). 
Solid lines are from the theory in the stable phase, continued as dashed lines into 
the unstable phases. Symbols are from simulations, circles, squares and diamonds are 
F 1,0,-1 respectively {N > 200, > 20 samples, > 10000 iterations). The insets 
show the case a = 0.2, F = — 1 for comparison. 



apparent discontinuity of (f) will become even more pronounced in the context of co- 
operation pressure, as discussed below. 



Statistical mechanics of a model eco-system 



15 



6. Effects of co-operation pressure 

We now turn to the effects of tlie co-operation pressure u on the behaviour of the model. 
We again limit the discussion to the case a = q = 1, and consider the system both with 
and without direct species interaction. 

6.1. No direct species interaction 

The phase behaviour of the system at w = is depicted as a function of the co-operation 
pressure in Fig. [HI For completeness we extend the discussion to positive and negative 
values of the co-operation pressure u, although only u > carries specific ecological 
meaning (plants however can grow without predation, which might be modelled by a 
positive self-interaction, corresponding to a negative co-operation pressure in the present 
model). Most interestingly the phase with optimal exploitation of resources is limited 
to the interval u = 0, a G [0, 0.27] on the u = axis. In particular, as mentioned 
above, any positive or negative amount of co-operation pressure removes the phase with 
H = 0. Furthermore as observed in Fig. [HI the eco-system is fully stable at all a for 
any positive co-operation pressure, even for infinitesimally small m > 0. For a > 0.27 
an unstable phase can only be found at m < Mc(a) < 0. Fig. [9] confirms that H > 
throughout this phase. As expected grows monotonically with u, the co-operation 
pressure u acts as a force driving the system into the interior of the simplex f[TOl) . For 
low or negative values of u on the other hand the fraction of surviving species is low. 
Our simulations seem to indicate that is continuous as m | for a > a^u = 0), 
but that a discontinuity may be present at lower values of a. This is similar to the 
behaviour of at low a a.s w I discussed above (see inset of Fig. [7j). As shown in 
the right panel of Fig. [9]0 attains values close to zero for a = 0.2 and u < 0, whereas 
the fraction of surviving species is clearly positive at positive u. While our simulations 
are potentially prone to finite-size effects, the data presented is consistent with a first 
order phase transition. Simulations furthermore indicate that might actually vanish 
at small enough a and negative co-operation pressure indicating the possible existence 
of a phase in which only a sub-extensive number of species survives. Such behaviour 
has previously been reported for the case of higher-order interaction in [TT] . Due to the 
limited relevance of negative co-operation pressure we have however not conducted a 
more detailed analysis of these observations, and can at this stage not fully confirm the 
existence of such a phase in this system of two-body interaction. 

6.2. With direct species interaction 

The phase structure of the model with co-operation pressure and direct species 
interaction is shown in Fig. [TOl Stable phases are found at large positive co-operation 
pressures and large relative numbers of resources, and either a reduction of u or a can 
induce instability. In line with earlier observations anti-symmetry in the direct species 
interactions tends to stabilise the eco-system, at full anti-correlation a stable fixed- 
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Figure 8. Phase diagram for tlie model witli co-operation pressure {w = 0, q = a = 1). 
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Figure 9. (Colour on-line) Effects of co-operation pressure: Efficiency of resource 
exploitation H (left) and fraction of surviving species cj) (right) versus co-operation 
pressure u(w = 0,(7 = (t = 1). Solid lines are from the theory in the stable phase, 
for a — 1, 1.5 continued as dashed lines into the unstable phases with finite integrated 
response. Symbols are from simulations (circles correspond to a = 0.2, squares to 
a — 1, diamonds to a = 1.5) with N = 300, run for 20000 iteration steps, averaged 
over 50 samples. Markers for a — 0.2 have been connected as a guide to the eye. 



point regime is found for any m > at any a, whereas unstable regimes can be found for 
r > —1 even at positive co-operation pressure. The left panel of Fig. [11] finally shows 
that H remains positive throughout all tested parameter ranges if w > 0. The right 
panel demonstrates that again (p is an increasing function of the co-operation pressure 
u. In contrast with the system at w = no discontinuities in the order parameters are 
observed, as the transition with diverging integrated response is absent. 
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Figure 10. (Colour on-line) Phase diagram for the model with species interaction 
{w — 1) and co-operation pressure. Resource variation is set to cr = 1. Curves show 
the onset of instability for F = 1, 0, —1 from top to bottom, with stable phases to the 
top-right, unstable ones to the lower left. 
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Figure 11. (Colour on-line) Effects of co-operation pressure: Efficiency of 
resource exploitation H (left) and fraction of surviving species (p (right) versus co- 
operation pressure u for model with direct species interaction (w — a = q = a = l). 
Solid lines are from the theory in the stable phase, continued as dashed lines into the 
unstable phases. Symbols are from simulations, circles correspond to F = 1, squares 
to F = 0, diamonds to F = — 1 (simulations are performed for N = 200, run for 10000 
iteration steps, averaged over 20 samples of the disorder). 



7. Effects of dilution 



Animals, of course, do not have the all-to-all interaction. We now turn to a discussion 
of the effects of diluting the interaction web between species. We restrict the discussion 
to the case without direct species interaction and without co-operation pressure, i.e. we 
consider only u = w = 0. We furthermore follow the philosophy of introducing dilution 
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in the context of neural networks [22] and in random replicator models [12] , and assume 
that only a fraction c G (0, 1] of interaction links between species is present. If an 
interaction between species i and j is present, then we take it to be determined by their 
respective use of resources, following [9]. More specifically, we write 

^^^'^-^E^'^'W + '^W (25) 



Xiit) N 



where 



N 

C, 



Qrw = Q{;-Ev^^^(^) (26) 



c 

3 



is the amount of resource ^ available to species i at time t. The coefficients Qj denote 
the dilution of the interactions and take values and 1, indicating a particular link to 
be absent or present, respectively. An absent link could for example correspond to a 
geographic separation between species. The Cij are here taken to be random, and we 
choose any Cij to be equal to one with probability c, and equal to with probability 
1 — c. Consequently we have 

c = (cii)c = (cii)c (27) 

for any pair i < j, where (. . denotes an average over realisations of the dilution. 
Note that {cfj)^ = c. Correlations in the interaction network are then introduced by 
the requirement that 

{cijCji)^ - = 7c(l - c) (28) 

with 7 G [0, 1]. 7 = 1 corresponds to an undirected symmetric network of interactions 
with Cij = Cji for all i < j. For 7 = Cij and Cji are uncorrected, and the links 
in the interaction web are hence directed ones. Ecologically realistic cases presumably 
correspond to 7 ~ 1, for completeness we extend the statistical mechanics analysis of the 
dilute model to general values 7 G [0, 1]. Finally, we note that following the conventions 
in the literature we write the number of resources as P = acN in this section, and 
that we take self-interactions to be present for all species, i.e. we have cu = 1 for all 
i = 1, . . . , N. We also note that we use Qq = a\fPC,^ and = along with the 



normalisation A^~^ 'Yli^ii^) ~ ^^^^ section 

The analysis of the dilute model is straightforward and can be performed along the 
lines of [29l [30], [12] . The effective process reads: 

x{t) = x{t) (^-a{l - c)x(t') - a ^ dt' [c{l - G)-^ + 7(1 - c)G] (t, t')x{t') - r]{t) + u{t) 

§ The modification to the statistics of the {<3o,Cf} is necessary to guarantee a well defined 
thermodynamic limit. While in the fully connected model all terms of order higher than A'^'^ drop 
out in the dynamical action due to the overall normalisation of species concentrations, this is would no 
longer be the case in the dilute model. If the statistics of the {Qo,Cf} were not modified, N different 
normalisation constraints would be required, due to different local interaction 'neighbourhoods' of 
species. The model specifications used in this section make sure that such terms do not appear. 
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Figure 12. (Colour on-line) Phase diagram for the dilute model (ct = 1). The 
curves show the transition lines below which fixed points become unstable. 7 = 
f , 0.75, 0.5, 0.25, from top to bottom. No divergence of the susceptibility x is observed 
for c < 1. At c = 1 one reproduces the transition of the model of [3], adc = 1) sa 0.27. 
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Figure 13. (Colour on-line) Effects of dilution: Fraction of surviving species (j) 
versus c for dilute model at a = 1.5 (cr =1). 7 = 0, 0.5, 1 from top to bottom. Solid 
lines are from theory in the stable phase, and have been continued as dashed lines 
into the unstable phase (where the theory can no longer be expected to be accurate). 
Symbols are from simulations {N = 300, run for 40000 time steps, averages over at 
least 10 samples are taken; for small values of c simulations may exhibit finite-size 
effects; also equilibration effects and sample to sample fluctuations cannot fully be 
excluded). The inset shows (P/Q versus c obtained from simulations for the same 
model parameters, and demonstrates the presence of a transition for 7 = 1 (circles). 
For 7 = 0, 0.5 d^/Q w in line with the theory which for those values of the symmetry 
parameter predicts the system to be stable for all c. 
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where 

{rjit)r]{t')) = a [c{l - G)-\aca^E + C)(I - G^)-^ + (1 - c)C] (t, t')(30) 

and with all other definitions as in the fully connected model. 

The resulting phase diagram is depicted in Fig. [121 As shown, correlated dilution 
(roughly 7 > 0.5) increases the numerical value of etc, and hence reduces the stable 
regime of the system compared with the fully connected model. At largely uncorrelated 
dilution 7 < 0.5 the location of the phase transition Oc shows only a weak dependence 
on the degree of dilution c. This behaviour is also refiected in Fig. [131 where we focus 
on the system at a = 1.5 and depict the fraction of surviving species as a function of 
the connectivity c at different values of the symmetry parameter 7. For 7 smaller than 
roughly one half, diluting the network of species does not seem to affect the stationary 
state significantly. The phase transition is absent, and the system always reaches a 
unique stable fixed point at this value of a, irrespectively of c. Uncorrelated dilution 
furthermore has only little effect on the diversity whereas as highly correlated interaction 
network can affect the ecosystem significantly, and reduces the number of survivors. 
This is in-line with our earlier observations on the effect of direct species interaction at 
different degrees of symmetry, see Fig. [Til 

8. Concluding remarks and outlook 

In summary we have used tools from disordered systems theory to study a stylised model 
of a simple eco-system, composed of a set of species competing for an amount of limited 
resources, and which at the same time are subject to direct inter-species competition. 

The dynamical system of corresponding replicator equations has been addressed 
by path integral techniques, allowing us in particular to study cases of asymmetric 
interaction between species (corresponding to prey-predator relations), where there is 
no Lyapunov function governing the dynamics, and where static approaches are hence 
inapplicable. 

We find that this simple model eco-system displays a rich spectrum of features, 
and interesting phase behaviour separating stable from unstable regimes. Our main 
findings can be summarised as follows: (i) in absence of direct species interaction and 
co-operation pressure the fully connected model displays a transition between a phase 
in which initial conditions are irrelevant and a non-ergodic phase. This transition is 
also marked by a change of the efficiency of resource exploitation. In the unstable 
phase resources are used optimally, while this is not the case in the stable phase, (ii) 
The introduction of either direct species interaction, co-operation pressure or dilution 
alters the type of transition observed, in particular the fully efficient phase is removed. 
One still finds a phase boundary separating a stable ergodic fixed point regime from 
a non-ergodic phase, (iii) For symmetric couplings the non-ergodic phase is marked 
by an exponential number of marginally stable fixed points, see also [HI [7] , and initial 
conditions determine which of these is reached in the long run leading to the observed 
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ergodicity breaking. At asymmetric or partially asymmetric coupling (induced by either 
direct interaction or dilution) no fixed point is reached in the unstable phase. Instead 
the trajectories of the system remain volatile and potentially chaotic, (iv) We observe a 
general tendency of increased stability when asymmetric (or anti-symmetric) interaction 
is introduced. This is the case both for direct species interaction and dilution. While 
the range of stability is then increased, no significant effects on the diversity of the 
eco-system are found, (v) The introduction of symmetric interaction or dilution can 
reduce the stability of the system significantly and at the same time also lead to a 
reduced diversity of its stationary population structure, (vi) In the absence of direct 
species interaction the effects of co-operation pressure can be drastic, and in particular 
the system is stable at any even infinitesimal amount of co-operation pressure. At small 
(relative) numbers of resources over species order parameters can display discontinuities 
as the co-operation pressure tends to zero. 

The model studied in the present paper and its phase behaviour are furthermore 
interesting from the statistical mechanics point of view. As detailed above the transition 
between a resource-efficient and an inefficient phase in the model without direct 
species interaction or co-operation pressure has a geometrical interpretation previously 
identified e.g. in the context of the dynamics of the Minority Game [261 EZj- This 
geometrical picture breaks down as soon as direct interaction or co-operation pressure 
are introduced, hence the absence of a transition at diverging integrated response and 
of the fully efficient phase. The present model may hence serve as a starting point for 
attempts to fully classify interacting agent models according to the presence or absence of 
phases with optimal resource exploitation. A close relation to the presence or otherwise 
of replica-symmetry breaking and to the geometry of the manifold of stationary states 
is here to be expected. 

Extensions of the present model might include adding a third or further trophic 
levels, temporally fluctuating resource availability (e.g. along the hues of [27]) or 
the introduction of further heterogeneity of the species. It is likely that this will 
alter the phase diagram, and might affect the stability or otherwise of the eco-system. 
Furthermore the computation of species abundance distributions (SAD), as introduced 
by Fisher et al |31] and by Preston [32] might be an interesting issue for future work. 
SAD have been measured and compared to log-normal and log-series distributions 
known in ecology for example in the model of [H] . The work of |33] demonstrates that 
random replicator models can yield SAD similar to left-skewed log-normal distributions. 
Given the presence of a phase transition in the model discussed here it would be 
particularly interesting to study finite systems near the transition, resulting in potential 
non-Gaussian features and fat-tailed abundance distributions (see [H [2] and references 
therein for similar critical fluctuations in Minority Game models near their phase 
transitions). In order to address the resulting topological structure and distribution 
of coupling strengths it might also be interesting to study the food web resulting from 
the present model in more detail. In particular, as seen above, some species die out 
asymptotically, inducing a reduced coupling matrix restricted to survivors only. While 
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species can not change their foraging strategies and no new hnks between species can 
created in our model, this extinction dynamics might lead to non-trivial, potentially 
correlated effective interaction strengths distributions among survivors at stationarity. 
In [3lj a dominance of prey-predator pairs in the set of surviving species has for example 
been identified in the context of replicator systems with quenched Gaussian interactions. 

Furthermore, it would be interesting to extend the analysis of the dilute model to 
more realistic finite-connectivity cases on complex networks [35J. So far we have only 
addressed dilute Erdos-Reyni type networks with an extensive number of connections 
per node. Complex networks with scale-free degree distribution [351 ES] and the small- 
world property [37] or other structures might here be of more biological relevance 
[m [151 HH]! in an approach to approximate dynamically evolved networks by static 
quenched ones. Extension to such sparse networks might require to study cases with 
only a finite number of interactions per species. This is challenging as the resulting 
effective dynamical theories do not close on the level of two-time order parameters. Still 
it would be interesting to examine the effects of network topology and degree sequence 
on the stability or otherwise of the model eco-system, as a first step potentially relying 
on numerical simulations or on replica approaches and the cavity-method [38] in order 
to study the statics of eco-systems with symmetric couplings. 
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